knitr::opts_chunk$set(echo = TRUE)

cmdstanr::set_cmdstan_path(path = "C:/Users/kueng/.cmdstan/cmdstan-2.35.0")


library(tidyverse)
library(R.utils)
library(wbCorr)
library(readxl)
library(kableExtra)
library(brms)
library(bayesplot)
library(beepr)
library(DHARMa)



source(file.path('Functions', 'ReportModels.R'))
source(file.path('Functions', 'PrettyTables.R'))
source(file.path('Functions', 'ReportMeasures.R'))
source(file.path('Functions', 'PrepareData.R'))
system("shutdown /a")
## [1] 1116
# Set options for analysis
use_mi = FALSE
shutdown = FALSE
report_ordinal = FALSE

options(
  dplyr.print_max = 100, 
  brms.backend = 'cmdstan',
  brms.file_refit = ifelse(use_mi, 'never', 'on_change'),
  brms.file_refit = 'on_change',
  #brms.file_refit = 'always',
  error = function() beepr::beep(sound = 5)
)



####################### Model parameters #######################

iterations = 5000
warmup = 2000
corstr = 'ar'
#corstr = 'cosy_couple'
#corstr = 'cosy_couple:user'


################################################################

if (corstr == 'ar') {
  autocor_str <- ~ ar(time = day, gr = coupleID:userID, p = 1)
  autocor_prior <- brms::set_prior("cauchy(0, 5)", class = "ar", lb = -1, ub = 1) 
  suffix = paste0('ar_', as.character(iterations))
} else if (corstr == 'cosy_couple') {
  autocor_str <- ~ cosy(gr = coupleID)
  autocor_prior <- brms::set_prior("cauchy(0, 5)", class = "cosy", lb = -1, ub = 1)
  suffix = paste0('cosy_couple_', as.character(iterations))
} else if (corstr == 'cosy_couple:user') {
  autocor_str <- ~ cosy(time = day, gr = coupleID:userID)
  autocor_prior <- brms::set_prior("cauchy(0, 5)", class = "cosy", lb = -1, ub = 1)
  suffix = paste0('cosy_couple:user_', as.character(iterations))
}
df <- openxlsx::read.xlsx(file.path('long.xlsx'))
df_original <- df

df_double <- prepare_data(df, recode_pushing = TRUE, use_mi = use_mi)[[1]]

Constructing scales Re-coding pusing reshaping data (4field) centering data within and between

Modelling

# For indistinguishable Dyads
model_rows_fixed <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'weartime_self_cb'
  )


model_rows_fixed_ordinal <- c(
  model_rows_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rows_fixed[2:length(model_rows_fixed)]
)

model_rows_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'ma[1]', 
  'cosy',
  'nu',
  'shape',
  'sderr',
  'sigma'
)

model_rows_random_ordinal <- c(model_rows_random,'disc')
# For indistinguishable Dyads
model_rownames_fixed <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'Daily received persuasion target -> target', 
    'Daily received persuasion target -> agent', 
    'Daily received pressure target -> target', 
    'Daily received pressure target -> agent', 
    'Daily received pushing target -> target', 
    'Daily received pushing target -> agent', 
    'Day', 
    'Daily weartime',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'Mean received persuasion target -> target',
    'Mean received persuasion target -> agent',
    'Mean received pressure target -> target',
    'Mean received pressure target -> agent',
    'Mean received pushing target -> target',
    'Mean received pushing target -> agent',
    'Mean weartime'
  )


model_rownames_fixed_ordinal <- c(
  model_rownames_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed[2:length(model_rownames_fixed)]
)

model_rownames_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(Daily received persuasion target -> target)', 
  'sd(Daily received persuasion target -> agent)', 
  'sd(Daily received pressure target -> target)', 
  'sd(Daily received pressure target -> agent)', 
  'sd(Daily received pushing target -> target)', 
  'sd(Daily received pushing target -> agent)', 
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]',
  'ma[1]', 
  'cosy',
  'nu',
  'shape',
  'sderr',
  'sigma'
)

model_rownames_random_ordinal <- c(model_rownames_random,'disc')
# For indistinguishable Dyads
model_rownames_fixed <- c(
    "Intercept", 
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily persuasion experienced", 
    "Daily persuasion utilized (partner's view)", # OR partner received
    "Daily pressure experienced", 
    "Daily pressure utilized (partner's view)", 
    "Daily pushing experienced", 
    "Daily pushing utilized (partner's view)", 
    "Day", 
    "Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean persuasion experienced", 
    "Mean persuasion utilized (partner's view)", 
    "Mean pressure experienced", 
    "Mean pressure utilized (partner's view)", 
    "Mean pushing experienced", 
    "Mean pushing utilized (partner's view)", 
    "Mean weartime"
  )


model_rownames_fixed_ordinal <- c(
  model_rownames_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed[2:length(model_rownames_fixed)]
)

model_rownames_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  "sd(Daily persuasion experienced)", 
  "sd(Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Daily pressure experienced)", 
  "sd(Daily pressure utilized (partner's view))", 
  "sd(Daily pushing experienced)", 
  "sd(Daily pushing utilized (partner's view))", 
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'ma[1]', 
  'cosy',
  'nu',
  'shape',
  'sderr',
  'sigma'
)

model_rownames_random_ordinal <- c(model_rownames_random,'disc')
rows_to_pack <- list(
  "Within-Person Effects" = c(2,9),
  "Between-Person Effects" = c(10,16),
  "Random Effects" = c(17, 23), 
  "Additional Parameters" = c(24,30)
  )


rows_to_pack_ordinal <- list(
  "Intercepts" = c(1,6),
  "Within-Person Effects" = c(2+5,9+5),
  "Between-Person Effects" = c(10+5,16+5),
  "Random Effects" = c(17+5, 23+5), 
  "Additional Parameters" = c(24+5,30+6)
  )

HURDLE MODELS

# For indistinguishable Dyads
model_rows_fixed_hu <- c(
    'Intercept', 
    'hu_Intercept',
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'weartime_self_cb',
    
    # HURDLE MODEL
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'hu_persuasion_self_cw', 
    'hu_persuasion_partner_cw', 
    'hu_pressure_self_cw', 
    'hu_pressure_partner_cw', 
    'hu_pushing_self_cw', 
    'hu_pushing_partner_cw', 
    'hu_day', 
    'hu_weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'hu_persuasion_self_cb',
    'hu_persuasion_partner_cb',
    'hu_pressure_self_cb',
    'hu_pressure_partner_cb',
    'hu_pushing_self_cb',
    'hu_pushing_partner_cb',
    'hu_weartime_self_cb'
  )

model_rows_random_hu <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(hu_Intercept)',
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # HURDLE
  'sd(hu_persuasion_self_cw)',
  'sd(hu_persuasion_partner_cw)',
  'sd(hu_pressure_self_cw)',
  'sd(hu_pressure_partner_cw)',
  'sd(hu_pushing_self_cw)',
  'sd(hu_pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'ma[1]', 
  'cosy',
  'nu',
  'shape',
  'sderr',
  'sigma'
)
# For indistinguishable Dyads
model_rownames_fixed_hu <- c(
    "Intercept", 
    "Hurdle Intercept",
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily persuasion experienced", 
    "Daily persuasion utilized (partner's view)", # OR partner received
    "Daily pressure experienced", 
    "Daily pressure utilized (partner's view)", 
    "Daily pushing experienced", 
    "Daily pushing utilized (partner's view)", 
    "Day", 
    "Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean persuasion experienced", 
    "Mean persuasion utilized (partner's view)", 
    "Mean pressure experienced", 
    "Mean pressure utilized (partner's view)", 
    "Mean pushing experienced", 
    "Mean pushing utilized (partner's view)", 
    "Mean weartime",
    
    # HURDLE
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Hu Daily persuasion experienced", 
    "Hu Daily persuasion utilized (partner's view)", # OR partner received
    "Hu Daily pressure experienced", 
    "Hu Daily pressure utilized (partner's view)", 
    "Hu Daily pushing experienced", 
    "Hu Daily pushing utilized (partner's view)", 
    "Hu Day", 
    "Hu Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Hu Mean persuasion experienced", 
    "Hu Mean persuasion utilized (partner's view)", 
    "Hu Mean pressure experienced", 
    "Hu Mean pressure utilized (partner's view)", 
    "Hu Mean pushing experienced", 
    "Hu Mean pushing utilized (partner's view)", 
    "Hu Mean weartime"
  )


model_rownames_random_hu <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(Hurdle Intercept)', 
  "sd(Daily persuasion experienced)", 
  "sd(Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Daily pressure experienced)", 
  "sd(Daily pressure utilized (partner's view))", 
  "sd(Daily pushing experienced)", 
  "sd(Daily pushing utilized (partner's view))", 
  
  # Hurdle
  "sd(Hu Daily persuasion experienced)", 
  "sd(Hu Daily persuasion utilized (partner's view))", # OR partner received
  "sd(Hu Daily pressure experienced)", 
  "sd(Hu Daily pressure utilized (partner's view))", 
  "sd(Hu Daily pushing experienced)", 
  "sd(Hu Daily pushing utilized (partner's view))", 
  # '-- CORRELATION STRUCTURE -- ', 
  'ar[1]', 
  'ma[1]', 
  'cosy',
  'nu',
  'shape',
  'sderr',
  'sigma'
)
rows_to_pack_hu <- list(
  "Conditional Within-Person Effects" = c(3,10),
  "Conditional Between-Person Effects" = c(11,17),
  
  "Hurdle Within-Person Effects" = c(18,25),
  "Hurdle Between-Person Effects" = c(26,32),
  
  "Random Effects" = c(33, 46), 
  "Additional Parameters" = c(47,53)
  )

Subjective MVPA

range(df_double$pa_sub, na.rm = T) 
## [1]   0 720
hist(df_double$pa_sub, breaks = 40) 

hist(log(df_double$pa_sub+00000000001), breaks = 40)

poisson

formula <- bf(
  pa_sub ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID),
  
  hu = ~ persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  , autocor = autocor_str
) 

prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 50)", class = "Intercept") # for non-zero PA
  , brms::set_prior("normal(6, 2.5)", class = "Intercept", dpar = 'hu') # hurdle part
  
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  
  #, brms::set_prior("normal(10, 10", class = "shape")
  
  #, brms::set_prior("cauchy(0, 2)", class = "sigma", lb = 0)
  , brms::set_prior("cauchy(0, 2)", class = "sderr", lb = 0)
  
  , autocor_prior
)

brms::validate_prior(
  prior1, 
  formula = formula, 
  data = df_double, 
  family = hurdle_poisson()
  )
## Warning: Rows containing NAs were excluded from the model.
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_sub <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  #family = brms::hurdle_lognormal(), 
  #family = brms::hurdle_negbinomial(), 
  family = brms::hurdle_poisson(),
  control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("pa_sub_hu_poisson", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
check_brms(pa_sub, log_pp_check = TRUE)
## 
## Divergences:
## 0 of 12000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 12000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 1672 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## 
## Computed from 12000 by 3736 log-likelihood matrix.
## 
##          Estimate     SE
## elpd_loo -32658.9  941.7
## p_loo     24314.6  853.7
## looic     65317.7 1883.5
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.4, 2.3]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     2064  55.2%   461     
##    (0.7, 1]   (bad)        98   2.6%   <NA>    
##    (1, Inf)   (very bad) 1574  42.1%   <NA>    
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(pa_sub, integer = TRUE)
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

## Object of Class DHARMa with simulated residuals based on 1000 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.3940511 1 0.02004498 0.308 0.5042459 0.925151 0.0009056866 0.4285656 0.6044282 0.5345751 0.9387669 0.1349345 0.5601726 1 0.4076557 0.9307523 0.3332225 0.3040894 0.09065378 0.998 ...
summarize_brms(
  pa_sub, 
  model_rows_fixed = model_rows_fixed_hu,
  model_rows_random = model_rows_random_hu,
  model_rownames_fixed = model_rownames_fixed_hu,
  model_rownames_random = model_rownames_random_hu,
  exponentiate = T,
  include_p_direction = T) %>%
  print_df(rows_to_pack = rows_to_pack_hu
)
## Warning in summarize_brms(pa_sub, model_rows_fixed = model_rows_fixed_hu, :
## Coefficients were exponentiated. Double check if this was intended.
exp(Est.) l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS p_direction
Intercept 48.82* 42.95 55.32 1.002 2631.97 4860.85 1
Hurdle Intercept 1.19 0.86 1.65 1.002 2447.79 5367.21 0.864
Conditional Within-Person Effects
Daily persuasion experienced 1.03 0.98 1.08 1.000 3873.02 6131.82 0.882
Daily persuasion utilized (partner’s view) 1.03 0.99 1.08 1.001 3448.17 5705.69 0.909
Daily pressure experienced 0.90* 0.82 0.99 1.000 3738.20 6032.04 0.98
Daily pressure utilized (partner’s view) 0.94 0.85 1.02 1.000 3797.25 5386.71 0.936
Daily pushing experienced 1.03 0.96 1.11 1.000 4416.62 6529.57 0.796
Daily pushing utilized (partner’s view) 0.99 0.93 1.06 1.001 3818.52 5967.44 0.56
Day 0.98 0.86 1.12 1.001 1900.30 3974.42 0.61
Daily weartime NA NA NA NA NA NA NA
Conditional Between-Person Effects
Mean persuasion experienced 1.03 0.75 1.42 1.001 2289.48 3858.10 0.577
Mean persuasion utilized (partner’s view) 1.02 0.74 1.39 1.001 2196.01 4213.69 0.539
Mean pressure experienced 1.09 0.75 1.57 1.001 2365.18 5282.43 0.677
Mean pressure utilized (partner’s view) 0.89 0.61 1.31 1.001 2776.35 5070.58 0.725
Mean pushing experienced 1.31 0.82 2.12 1.001 2771.56 4517.09 0.874
Mean pushing utilized (partner’s view) 1.31 0.81 2.14 1.001 2628.55 4641.10 0.862
Mean weartime NA NA NA NA NA NA NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced 0.65* 0.57 0.73 1.001 10583.80 8967.97 1
Hu Daily persuasion utilized (partner’s view) 0.75* 0.67 0.84 1.000 11285.52 9314.41 1
Hu Daily pressure experienced 1.23 0.89 1.72 1.000 11630.87 9358.50 0.904
Hu Daily pressure utilized (partner’s view) 0.66* 0.43 0.95 1.000 9797.07 6380.89 0.987
Hu Daily pushing experienced 0.58* 0.40 0.78 1.000 9744.20 8381.38 1
Hu Daily pushing utilized (partner’s view) 0.54* 0.41 0.68 1.001 10961.44 8955.36 1
Hu Day 1.09 0.84 1.42 1.001 23427.12 8694.42 0.75
Hu Daily weartime NA NA NA NA NA NA NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced 0.82 0.35 1.95 1.001 2043.52 4080.81 0.676
Hu Mean persuasion utilized (partner’s view) 0.83 0.36 1.95 1.001 2064.02 4271.06 0.668
Hu Mean pressure experienced 3.49* 1.40 8.69 1.000 3324.64 6122.70 0.996
Hu Mean pressure utilized (partner’s view) 1.84 0.74 4.54 1.001 3237.84 5819.30 0.904
Hu Mean pushing experienced 0.32 0.09 1.10 1.001 3339.43 5588.58 0.965
Hu Mean pushing utilized (partner’s view) 0.31 0.09 1.09 1.001 3358.06 5136.34 0.966
Hu Mean weartime NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.31 0.24 0.41 1.00 2997.87 5058.08 NA
sd(Hurdle Intercept) 0.90 0.70 1.17 1.00 3564.02 5744.88 NA
sd(Daily persuasion experienced) 0.11 0.07 0.16 1.00 2668.43 4480.13 NA
sd(Daily persuasion utilized (partner’s view)) 0.08 0.04 0.13 1.00 2280.60 3048.95 NA
sd(Daily pressure experienced) 0.06 0.00 0.20 1.00 2887.46 4982.02 NA
sd(Daily pressure utilized (partner’s view)) 0.06 0.00 0.18 1.00 2709.88 5297.01 NA
sd(Daily pushing experienced) 0.12 0.06 0.20 1.00 2030.79 2407.90 NA
sd(Daily pushing utilized (partner’s view)) 0.11 0.04 0.18 1.00 1884.95 2419.13 NA
sd(Hu Daily persuasion experienced) 0.18 0.02 0.34 1.00 3860.73 3423.70 NA
sd(Hu Daily persuasion utilized (partner’s view)) 0.17 0.02 0.33 1.00 3691.94 3266.04 NA
sd(Hu Daily pressure experienced) 0.30 0.01 0.85 1.00 4749.58 6186.83 NA
sd(Hu Daily pressure utilized (partner’s view)) 0.35 0.01 1.02 1.00 4522.11 5369.01 NA
sd(Hu Daily pushing experienced) 0.65 0.33 1.08 1.00 6169.71 8308.81 NA
sd(Hu Daily pushing utilized (partner’s view)) 0.32 0.05 0.65 1.00 5399.56 4536.50 NA
Additional Parameters
ar[1] 0.35 0.27 0.42 1.01 823.84 1687.17 NA
ma[1] NA NA NA NA NA NA NA
cosy NA NA NA NA NA NA NA
nu NA NA NA NA NA NA NA
shape NA NA NA NA NA NA NA
sderr 0.62 0.60 0.65 1.00 2228.33 4497.83 NA
sigma NA NA NA NA NA NA NA

Device Based MVPA

range(df_double$pa_obj, na.rm = T) 
## [1]   5.75 971.25
hist(df_double$pa_obj, breaks = 50)

df_double$pa_obj_log <- log(df_double$pa_obj)

hist(df_double$pa_obj_log, breaks = 50)

lognormal

formula <- bf(
  pa_obj ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  , autocor = autocor_str
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 50)", class = "Intercept") # for non-zero PA

  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  
  , brms::set_prior("cauchy(0, 2)", class = "sigma", lb = 0)
  , brms::set_prior("cauchy(0, 2)", class = "sderr", lb = 0)
  
  , autocor_prior
)

brms::validate_prior(
  prior1, 
  formula = formula, 
  data = df_double, 
  family = lognormal()
  )
## Warning: Rows containing NAs were excluded from the model.
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_obj_log <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = lognormal(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("pa_obj_lognormal_", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
check_brms(pa_obj_log, log_pp_check = TRUE)
## 
## Divergences:
## 12 of 12000 iterations ended with a divergence (0.1%).
## Try increasing 'adapt_delta' to remove the divergences.
## 
## Tree depth:
## 2114 of 12000 iterations saturated the maximum tree depth of 10 (17.6166666666667%).
## Try increasing 'max_treedepth' to avoid saturation.
## 
## Energy:
## E-BFMI indicated possible pathological behavior:
##   Chain 1: E-BFMI = 0.007
##   Chain 2: E-BFMI = 0.007
##   Chain 3: E-BFMI = 0.012
##   Chain 4: E-BFMI = 0.005
## E-BFMI below 0.2 indicates you may need to reparameterize your model.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 2959 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## 
## Computed from 12000 by 3337 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -20147.3 105.0
## p_loo      3902.6  80.7
## looic     40294.6 210.1
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.0, 1.4]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)      378  11.3%   3       
##    (0.7, 1]   (bad)      2424  72.6%   <NA>    
##    (1, Inf)   (very bad)  535  16.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(pa_obj_log, integer = TRUE)
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

## Object of Class DHARMa with simulated residuals based on 1000 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.177 0.826 0.004 0.024 0.082 0.825 0.162 0.251 0.39 0.201 0.67 0.262 0.167 0.285 0.048 0.552 0.375 0.318 0.642 0.706 ...
summarize_brms(
  pa_obj_log, 
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T,
  include_p_direction = T) %>%
  print_df(rows_to_pack = rows_to_pack)
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be
## careful when analysing the results! We recommend running more iterations and/or
## setting stronger priors.
## Warning: There were 12 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning in summarize_brms(pa_obj_log, model_rows_fixed = model_rows_fixed, :
## Coefficients were exponentiated. Double check if this was intended.
exp(Est.) l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS p_direction
Intercept 117.78* 105.84 130.68 1.001 1601.54 2924.15 1
Within-Person Effects
Daily persuasion experienced 1.03 1.00 1.06 1.001 3821.19 5697.74 0.975
Daily persuasion utilized (partner’s view) 1.02 0.99 1.05 1.001 4680.21 6890.35 0.888
Daily pressure experienced 0.95 0.89 1.01 1.000 4269.96 6167.77 0.949
Daily pressure utilized (partner’s view) 0.98 0.92 1.04 1.001 4783.85 6711.79 0.723
Daily pushing experienced 1.03 0.98 1.08 1.002 4357.67 5639.41 0.89
Daily pushing utilized (partner’s view) 1.03 0.99 1.06 1.000 4484.91 5950.21 0.903
Day 0.97 0.89 1.05 1.001 3479.74 5345.43 0.796
Daily weartime 1.00* 1.00 1.00 1.000 5566.97 8054.65 1
Between-Person Effects
Mean persuasion experienced 1.08 0.81 1.45 1.004 1294.98 2335.94 0.714
Mean persuasion utilized (partner’s view) 0.96 0.72 1.29 1.004 1269.93 1823.69 0.603
Mean pressure experienced 0.98 0.72 1.32 1.002 1745.67 3666.00 0.561
Mean pressure utilized (partner’s view) 0.98 0.73 1.32 1.004 1420.53 2319.05 0.548
Mean pushing experienced 1.00 0.66 1.51 1.003 1941.87 4136.41 0.509
Mean pushing utilized (partner’s view) 1.27 0.84 1.92 1.002 2077.74 4185.41 0.876
Mean weartime 1.00 1.00 1.00 1.001 4632.74 6314.23 0.878
Random Effects
sd(Intercept) 0.30 0.23 0.39 1.00 1842.51 3748.05 NA
sd(Daily persuasion experienced) 0.05 0.03 0.09 1.00 2705.63 3893.37 NA
sd(Daily persuasion utilized (partner’s view)) 0.05 0.02 0.08 1.00 2276.55 2313.82 NA
sd(Daily pressure experienced) 0.04 0.00 0.13 1.00 2192.13 3758.21 NA
sd(Daily pressure utilized (partner’s view)) 0.04 0.00 0.11 1.00 3285.19 5044.16 NA
sd(Daily pushing experienced) 0.06 0.00 0.14 1.01 899.66 2385.76 NA
sd(Daily pushing utilized (partner’s view)) 0.03 0.00 0.08 1.00 2257.27 4488.51 NA
Additional Parameters
ar[1] 0.39 0.26 0.74 1.22 16.14 113.59 NA
ma[1] NA NA NA NA NA NA NA
cosy NA NA NA NA NA NA NA
nu NA NA NA NA NA NA NA
shape NA NA NA NA NA NA NA
sderr 0.45 0.32 0.53 1.25 14.86 82.72 NA
sigma 0.30 0.18 0.42 1.25 14.80 63.26 NA

log gaussian

formula <- bf(
  log(pa_obj) ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  , autocor = autocor_str
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 50)", class = "Intercept") # for non-zero PA

  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  
  , brms::set_prior("cauchy(0, 2)", class = "sigma", lb = 0)
  #, brms::set_prior("cauchy(0, 2)", class = "sderr", lb = 0)
  
  , autocor_prior
)

brms::validate_prior(
  prior1, 
  formula = formula, 
  data = df_double, 
  family = gaussian()
  )
## Warning: Rows containing NAs were excluded from the model.
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_obj_log <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian_", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
check_brms(pa_obj_log, log_pp_check = TRUE)
## 
## Divergences:
## 0 of 12000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 12000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## 
## Computed from 12000 by 3337 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -2813.8  56.5
## p_loo        91.6   4.6
## looic      5627.6 113.0
## ------
## MCSE of elpd_loo is 0.1.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.3, 1.7]).
## 
## All Pareto k estimates are good (k < 0.7).
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(pa_obj_log, integer = TRUE)
## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

## DHARMa:testOutliers with type = binomial may have inflated Type I error rates for integer-valued distributions. To get a more exact result, it is recommended to re-run testOutliers with type = 'bootstrap'. See ?testOutliers for details

## Object of Class DHARMa with simulated residuals based on 1000 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.379 0.677 0.011 0.118 0.106 0.568 0.206 0.574 0.267 0.296 0.761 0.514 0.395 0.524 0.086 0.381 0.398 0.255 0.395 0.729 ...
summarize_brms(
  pa_obj_log, 
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T,
  include_p_direction = T) %>%
  print_df(rows_to_pack = rows_to_pack)
## Warning in summarize_brms(pa_obj_log, model_rows_fixed = model_rows_fixed, :
## Coefficients were exponentiated. Double check if this was intended.
exp(Est.) l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS p_direction
Intercept 117.92* 105.70 131.58 1.000 3121.13 5697.81 1
Within-Person Effects
Daily persuasion experienced 1.03 1.00 1.06 1.000 8124.19 8894.78 0.958
Daily persuasion utilized (partner’s view) 1.02 0.99 1.05 1.000 10231.49 8773.49 0.878
Daily pressure experienced 0.95 0.89 1.01 1.001 12898.90 8806.28 0.953
Daily pressure utilized (partner’s view) 0.98 0.92 1.04 1.000 13958.09 8679.13 0.731
Daily pushing experienced 1.03 0.98 1.07 1.000 9673.48 7604.96 0.881
Daily pushing utilized (partner’s view) 1.03 0.99 1.07 1.000 13420.04 8709.55 0.943
Day 0.96 0.88 1.05 1.000 16584.01 8216.44 0.817
Daily weartime 1.00* 1.00 1.00 1.001 11954.84 7895.70 1
Between-Person Effects
Mean persuasion experienced 1.09 0.82 1.45 1.001 2669.09 4389.54 0.72
Mean persuasion utilized (partner’s view) 0.96 0.71 1.29 1.002 2642.26 4329.30 0.591
Mean pressure experienced 0.97 0.71 1.34 1.001 4130.93 6947.20 0.584
Mean pressure utilized (partner’s view) 0.98 0.72 1.33 1.001 4031.90 6870.09 0.56
Mean pushing experienced 1.01 0.66 1.54 1.001 4026.01 5657.99 0.518
Mean pushing utilized (partner’s view) 1.29 0.85 1.96 1.001 4011.26 6538.24 0.886
Mean weartime 1.00 1.00 1.00 1.001 11714.34 9448.27 0.836
Random Effects
sd(Intercept) 0.30 0.23 0.39 1.00 3616.86 6285.29 NA
sd(Daily persuasion experienced) 0.05 0.03 0.08 1.00 6572.95 6506.41 NA
sd(Daily persuasion utilized (partner’s view)) 0.05 0.02 0.08 1.00 5974.21 5245.51 NA
sd(Daily pressure experienced) 0.04 0.00 0.13 1.00 5894.74 5871.09 NA
sd(Daily pressure utilized (partner’s view)) 0.04 0.00 0.11 1.00 5855.56 5363.44 NA
sd(Daily pushing experienced) 0.06 0.00 0.14 1.00 2348.35 4924.03 NA
sd(Daily pushing utilized (partner’s view)) 0.03 0.00 0.08 1.00 4996.80 6091.79 NA
Additional Parameters
ar[1] 0.29 0.26 0.33 1.00 17999.27 8503.02 NA
ma[1] NA NA NA NA NA NA NA
cosy NA NA NA NA NA NA NA
nu NA NA NA NA NA NA NA
shape NA NA NA NA NA NA NA
sderr NA NA NA NA NA NA NA
sigma 0.55 0.54 0.57 1.00 16959.11 9077.14 NA

Affect

range(df_double$aff, na.rm = T) 
## [1] 0 5
hist(df_double$aff, breaks = 15)

Gaussian

formula <- bf(
  aff ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  , autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b")
  ,brms::set_prior("normal(0, 20)", class = "Intercept", lb=1, ub=6) # range of the outcome scale
  
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  
  , brms::set_prior("cauchy(0, 2)", class = "sigma", lb = 0)
  #, brms::set_prior("cauchy(0, 2)", class = "sderr", lb = 0)
  
  , autocor_prior
)

brms::validate_prior(
  prior1, 
  formula = formula, 
  data = df_double, 
  family = gaussian()
  )
## Warning: Rows containing NAs were excluded from the model.
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

mood <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("mood_gauss_", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
check_brms(mood, log_pp_check = FALSE)
## 
## Divergences:
## 0 of 12000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 12000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 1 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## 
## Computed from 12000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo  -4816.5  63.5
## p_loo        89.9   4.6
## looic      9633.1 127.1
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.5, 1.8]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     3735  100.0%  1040    
##    (0.7, 1]   (bad)         1    0.0%  <NA>    
##    (1, Inf)   (very bad)    0    0.0%  <NA>    
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(mood, integer = FALSE)

## Object of Class DHARMa with simulated residuals based on 1000 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.504 0.804 0.261 0.971 0.295 0.576 0.358 0.331 0.373 0.468 0.607 0.056 0.657 0.828 0.443 0.496 0.785 0.76 0.629 0.485 ...
summarize_brms(
  mood, 
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = F,
  include_p_direction = T) %>%
  print_df(rows_to_pack = rows_to_pack)
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS p_direction
Intercept 3.73* 3.52 3.94 1.001 4370.23 6385.35 1
Within-Person Effects
Daily persuasion experienced 0.00 -0.03 0.04 1.000 17939.10 9789.59 0.536
Daily persuasion utilized (partner’s view) 0.02 -0.02 0.06 1.000 16083.36 9736.51 0.792
Daily pressure experienced -0.05 -0.16 0.05 1.001 13407.68 8929.88 0.839
Daily pressure utilized (partner’s view) -0.04 -0.18 0.09 1.000 12032.49 8988.62 0.712
Daily pushing experienced 0.02 -0.04 0.09 1.001 14740.42 10175.51 0.761
Daily pushing utilized (partner’s view) 0.06* 0.01 0.12 1.000 15726.84 9585.61 0.983
Day 0.22* 0.06 0.38 1.000 23524.10 8669.26 0.996
Daily weartime NA NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 0.31 -0.26 0.89 1.001 3497.64 5376.58 0.859
Mean persuasion utilized (partner’s view) 0.23 -0.35 0.78 1.001 3509.40 5215.53 0.789
Mean pressure experienced -0.29 -0.88 0.30 1.001 4509.37 6827.12 0.833
Mean pressure utilized (partner’s view) -0.31 -0.89 0.28 1.001 4639.37 6766.96 0.849
Mean pushing experienced 0.25 -0.55 1.05 1.001 5563.03 7538.13 0.73
Mean pushing utilized (partner’s view) 0.38 -0.41 1.19 1.000 5684.30 7991.92 0.832
Mean weartime NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.59 0.46 0.77 1.00 4200.42 6902.14 NA
sd(Daily persuasion experienced) 0.03 0.00 0.07 1.00 6399.64 7211.47 NA
sd(Daily persuasion utilized (partner’s view)) 0.06 0.01 0.11 1.00 3965.10 4144.64 NA
sd(Daily pressure experienced) 0.11 0.01 0.28 1.00 4670.61 6951.38 NA
sd(Daily pressure utilized (partner’s view)) 0.19 0.03 0.38 1.00 4325.92 4766.79 NA
sd(Daily pushing experienced) 0.09 0.02 0.17 1.00 5567.01 4535.34 NA
sd(Daily pushing utilized (partner’s view)) 0.05 0.00 0.13 1.00 5899.70 6469.91 NA
Additional Parameters
ar[1] 0.45 0.42 0.48 1.00 21354.19 8810.71 NA
ma[1] NA NA NA NA NA NA NA
cosy NA NA NA NA NA NA NA
nu NA NA NA NA NA NA NA
shape NA NA NA NA NA NA NA
sderr NA NA NA NA NA NA NA
sigma 0.87 0.85 0.89 1.00 21666.82 8355.30 NA

Skewnormal

formula <- bf(
  aff ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  , autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b")
  ,brms::set_prior("normal(0, 20)", class = "Intercept", lb=1, ub=6) # range of the outcome scale
  
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  
  , brms::set_prior("cauchy(0, 2)", class = "sigma", lb = 0)
  , brms::set_prior("cauchy(0, 2)", class = "sderr", lb = 0)
  
  , autocor_prior
)

brms::validate_prior(
  prior1, 
  formula = formula, 
  data = df_double, 
  family = brms::skew_normal()
  )
## Warning: Rows containing NAs were excluded from the model.
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

mood <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::skew_normal(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("mood_skew_", suffix))
)
## Warning: Rows containing NAs were excluded from the model.
check_brms(mood, log_pp_check = FALSE)
## 
## Divergences:
## 0 of 12000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 12000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 1495 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## 
## Computed from 12000 by 3736 log-likelihood matrix.
## 
##          Estimate    SE
## elpd_loo -10406.6 341.7
## p_loo      6135.2 332.1
## looic     20813.2 683.4
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.2, 1.5]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     2241  60.0%   584     
##    (0.7, 1]   (bad)       145   3.9%   <NA>    
##    (1, Inf)   (very bad) 1350  36.1%   <NA>    
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(mood, integer = FALSE)

## Object of Class DHARMa with simulated residuals based on 1000 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.204 0.656 0.269 0.602 0.329 0.452 0.435 0.22 0.272 0.402 0.569 0.077 0.407 0.714 0.394 0.387 0.926 0.737 0.601 0.289 ...
summarize_brms(
  mood, 
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = F,
  include_p_direction = T) %>%
  print_df(rows_to_pack = rows_to_pack)
b l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS p_direction
Intercept 3.73* 3.58 3.88 1.002 1373.37 2791.19 1
Within-Person Effects
Daily persuasion experienced -0.01 -0.03 0.02 1.000 8600.23 8725.27 0.693
Daily persuasion utilized (partner’s view) 0.01 -0.01 0.04 1.000 8312.57 8754.48 0.811
Daily pressure experienced -0.02 -0.09 0.05 1.000 7981.67 8792.42 0.751
Daily pressure utilized (partner’s view) -0.02 -0.10 0.05 1.000 6427.79 6194.41 0.753
Daily pushing experienced 0.03 -0.01 0.07 1.000 7934.70 7875.14 0.918
Daily pushing utilized (partner’s view) 0.04 -0.01 0.08 1.000 8388.32 8465.48 0.947
Day 0.08 0.00 0.16 1.000 5062.03 7646.57 0.971
Daily weartime NA NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 0.23 -0.17 0.63 1.001 1283.06 2482.71 0.872
Mean persuasion utilized (partner’s view) 0.06 -0.34 0.46 1.001 1294.51 2610.65 0.619
Mean pressure experienced -0.13 -0.52 0.26 1.002 1430.26 3391.52 0.746
Mean pressure utilized (partner’s view) -0.20 -0.58 0.20 1.001 1512.29 3250.52 0.845
Mean pushing experienced 0.02 -0.54 0.59 1.001 1983.70 3529.67 0.537
Mean pushing utilized (partner’s view) 0.37 -0.19 0.94 1.001 1967.50 3588.70 0.904
Mean weartime NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.42 0.33 0.55 1.00 2293.20 3885.79 NA
sd(Daily persuasion experienced) 0.02 0.00 0.05 1.00 3349.92 4252.26 NA
sd(Daily persuasion utilized (partner’s view)) 0.02 0.00 0.05 1.00 3283.81 4707.76 NA
sd(Daily pressure experienced) 0.04 0.00 0.14 1.00 5898.15 6001.61 NA
sd(Daily pressure utilized (partner’s view)) 0.06 0.00 0.17 1.00 3880.39 5131.95 NA
sd(Daily pushing experienced) 0.03 0.00 0.07 1.00 4684.33 5777.64 NA
sd(Daily pushing utilized (partner’s view)) 0.04 0.00 0.10 1.00 3374.46 3951.91 NA
Additional Parameters
ar[1] 0.90 0.74 1.00 1.00 4416.46 4790.19 NA
ma[1] NA NA NA NA NA NA NA
cosy NA NA NA NA NA NA NA
nu NA NA NA NA NA NA NA
shape NA NA NA NA NA NA NA
sderr 0.20 0.15 0.25 1.00 759.51 1614.13 NA
sigma 0.86 0.82 0.90 1.00 860.04 1834.85 NA

Reactance

Gaussian

range(df_double$reactance, na.rm = T) 
## [1] 0 5
hist(df_double$reactance, breaks = 6)

hist(log(df_double$reactance + 0.0000001))

formula <- bf(
  reactance ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID), 
  
  hu = ~ persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  , autocor = autocor_str
)



prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b")
  , brms::set_prior("normal(0, 20)", class = "Intercept") # for non-zero PA
  , brms::set_prior("normal(20, 10)", class = "Intercept", dpar = 'hu') # hurdle part
  
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  
  , brms::set_prior("cauchy(0, 2)", class = "sigma", lb = 0)
  , brms::set_prior("cauchy(0, 2)", class = "sderr", lb = 0)
  
  , autocor_prior
)


brms::validate_prior(
  prior1, 
  formula = formula, 
  data = df_double, 
  family = hurdle_lognormal()
  )
## Warning: Rows containing NAs were excluded from the model.
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::hurdle_lognormal(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("reactance_hu_", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
check_brms(reactance, log_pp_check = T)
## 
## Divergences:
## 330 of 12000 iterations ended with a divergence (2.75%).
## Try increasing 'adapt_delta' to remove the divergences.
## 
## Tree depth:
## 0 of 12000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 171 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## 
## Computed from 12000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -751.8 35.9
## p_loo       247.6 16.7
## looic      1503.6 71.8
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.0, 1.3]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     585   77.4%   15      
##    (0.7, 1]   (bad)      151   20.0%   <NA>    
##    (1, Inf)   (very bad)  20    2.6%   <NA>    
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(reactance, integer = FALSE)

## Object of Class DHARMa with simulated residuals based on 1000 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.5810832 0.4374127 0.841 0.08129721 0.6129632 0.941 0.6289292 0.1092773 0.938 0.03320617 0.07669713 0.1846535 0.936 0.4611895 0.931 0.7176145 0.3412033 0.7992349 0.09020157 0.4778683 ...
summarize_brms(
  reactance, 
  model_rows_fixed = model_rows_fixed_hu,
  model_rows_random = model_rows_random_hu,
  model_rownames_fixed = model_rownames_fixed_hu,
  model_rownames_random = model_rownames_random_hu,
  exponentiate = T,
  include_p_direction = T) %>%
  print_df(rows_to_pack = rows_to_pack_hu)
## Warning: There were 330 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning in summarize_brms(reactance, model_rows_fixed = model_rows_fixed_hu, :
## Coefficients were exponentiated. Double check if this was intended.
exp(Est.) l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS p_direction
Intercept 1.80* 1.49 2.15 1.003 1110.33 314.21 1
Hurdle Intercept 4.59* 2.34 9.06 1.005 1116.76 260.40 1
Conditional Within-Person Effects
Daily persuasion experienced 0.99 0.92 1.07 1.003 5854.63 6400.86 0.598
Daily persuasion utilized (partner’s view) 1.00 0.93 1.06 1.001 9909.08 9180.09 0.537
Daily pressure experienced 1.06 0.93 1.18 1.001 5299.52 5745.59 0.852
Daily pressure utilized (partner’s view) 1.05 0.87 1.29 1.001 7738.57 6558.13 0.705
Daily pushing experienced 1.01 0.95 1.08 1.001 4837.74 4729.72 0.616
Daily pushing utilized (partner’s view) 1.01 0.90 1.15 1.000 5688.36 6664.54 0.543
Day 1.00 0.78 1.28 1.001 5590.68 5465.10 0.507
Daily weartime NA NA NA NA NA NA NA
Conditional Between-Person Effects
Mean persuasion experienced 0.94 0.65 1.36 1.002 1943.80 760.30 0.64
Mean persuasion utilized (partner’s view) 1.22 0.82 1.85 1.001 3824.99 6081.63 0.838
Mean pressure experienced 1.21 0.87 1.72 1.001 2625.99 3181.89 0.867
Mean pressure utilized (partner’s view) 0.86 0.60 1.23 1.002 3365.88 8354.95 0.784
Mean pushing experienced 0.85 0.52 1.36 1.001 7080.59 7313.42 0.756
Mean pushing utilized (partner’s view) 1.00 0.47 2.04 1.001 6533.12 7762.23 0.501
Mean weartime NA NA NA NA NA NA NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced 1.19 0.98 1.47 1.001 12763.62 8749.26 0.961
Hu Daily persuasion utilized (partner’s view) 0.96 0.71 1.31 1.000 7107.57 8188.27 0.625
Hu Daily pressure experienced 0.53 0.26 1.06 1.002 2641.96 2713.86 0.965
Hu Daily pressure utilized (partner’s view) 0.75 0.25 2.25 1.003 2617.19 3019.63 0.738
Hu Daily pushing experienced 0.80 0.63 1.01 1.002 12344.22 9031.91 0.968
Hu Daily pushing utilized (partner’s view) 1.14 0.78 1.75 1.002 10356.58 6278.80 0.751
Hu Day 0.61 0.28 1.30 1.001 17742.68 8649.30 0.898
Hu Daily weartime NA NA NA NA NA NA NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced 0.63 0.15 2.32 1.001 3139.60 1532.65 0.748
Hu Mean persuasion utilized (partner’s view) 0.59 0.13 2.53 1.002 5231.46 2460.22 0.757
Hu Mean pressure experienced 0.05* 0.00 0.50 1.001 7284.83 7847.82 0.995
Hu Mean pressure utilized (partner’s view) 0.84 0.08 10.66 1.002 2473.08 1860.09 0.559
Hu Mean pushing experienced 1.75 0.21 17.49 1.000 5968.83 7326.48 0.682
Hu Mean pushing utilized (partner’s view) 25.69* 2.54 336.83 1.001 6176.19 5147.43 0.997
Hu Mean weartime NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.17 0.01 0.36 1.00 1423.67 980.63 NA
sd(Hurdle Intercept) 1.12 0.67 1.73 1.00 2715.38 1741.14 NA
sd(Daily persuasion experienced) 0.11 0.04 0.20 1.00 2264.53 3904.62 NA
sd(Daily persuasion utilized (partner’s view)) 0.06 0.00 0.16 1.00 3632.58 5398.10 NA
sd(Daily pressure experienced) 0.12 0.01 0.28 1.00 4027.14 4485.55 NA
sd(Daily pressure utilized (partner’s view)) 0.15 0.01 0.47 1.00 2816.95 4334.15 NA
sd(Daily pushing experienced) 0.07 0.00 0.19 1.00 2239.61 5326.11 NA
sd(Daily pushing utilized (partner’s view)) 0.12 0.01 0.29 1.00 2702.50 6111.61 NA
sd(Hu Daily persuasion experienced) 0.22 0.01 0.51 1.00 2163.81 4638.10 NA
sd(Hu Daily persuasion utilized (partner’s view)) 0.49 0.08 0.99 1.00 2716.33 3237.37 NA
sd(Hu Daily pressure experienced) 1.11 0.13 2.45 1.00 2696.21 2624.71 NA
sd(Hu Daily pressure utilized (partner’s view)) 1.10 0.05 3.22 1.00 4013.21 6004.97 NA
sd(Hu Daily pushing experienced) 0.24 0.01 0.59 1.00 4288.87 5981.99 NA
sd(Hu Daily pushing utilized (partner’s view)) 0.32 0.01 0.99 1.00 3278.34 5905.15 NA
Additional Parameters
ar[1] 0.19 -0.85 0.94 1.00 1534.22 3457.71 NA
ma[1] NA NA NA NA NA NA NA
cosy NA NA NA NA NA NA NA
nu NA NA NA NA NA NA NA
shape NA NA NA NA NA NA NA
sderr 0.19 0.01 0.38 1.01 255.37 247.93 NA
sigma 0.33 0.18 0.44 1.01 248.57 182.60 NA 
hypothesis(reactance, "pressure_self_cw > pushing_self_cw")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (pressure_self_cw... > 0     0.05      0.07    -0.07     0.16       3.51
##   Post.Prob Star
## 1      0.78     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.
hypothesis(reactance, "hu_pressure_self_cw > hu_pushing_self_cw")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (hu_pressure_self... > 0    -0.41      0.39    -1.03     0.21       0.15
##   Post.Prob Star
## 1      0.13     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Ordinal

df_double$reactance_ordinal <- factor(df_double$reactance,
                                      levels = 0:5, 
                                      ordered = TRUE)

formula <- bf(
  reactance_ordinal ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  , autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  #, brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5) # range of the outcome scale
  
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  
  #, brms::set_prior("cauchy(0, 2)", class = "sigma", lb = 0)
  , brms::set_prior("cauchy(0, 2)", class = "sderr", lb = 0)
  
  , autocor_prior
)


brms::validate_prior(
  prior1, 
  formula = formula, 
  data = df_double, 
  family = cumulative() # HURDLE_CUMULATIVE
  )
## Warning: Rows containing NAs were excluded from the model.
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance_ordinal <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::cumulative(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("reactance_ordinal_noar_", suffix))
  , file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
check_brms(reactance_ordinal)
## 
## Divergences:
## 0 of 12000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 12000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 5 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## 
## Computed from 12000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -683.4 32.1
## p_loo        92.7  6.2
## looic      1366.9 64.2
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.3, 1.4]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)     751   99.3%   187     
##    (0.7, 1]   (bad)        5    0.7%   <NA>    
##    (1, Inf)   (very bad)   0    0.0%   <NA>    
## See help('pareto-k-diagnostic') for details.
summarize_brms(
  reactance_ordinal, 
  model_rows_fixed = model_rows_fixed_ordinal,
  model_rows_random = model_rows_random_ordinal,
  model_rownames_fixed = model_rownames_fixed_ordinal,
  model_rownames_random = model_rownames_random_ordinal,
  exponentiate = T,
  include_p_direction = T) %>%
  print_df(rows_to_pack = rows_to_pack_ordinal)
OR l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS p_direction
Intercepts
Intercept NA NA NA NA NA NA NA
Intercept[1] 4.03* 2.37 7.11 1.000 5771.42 6583.42 1
Intercept[2] 8.94* 5.13 16.70 1.001 4147.38 5353.11 1
Intercept[3] 25.53* 13.66 52.08 1.001 3636.11 3983.04 1
Intercept[4] 115.36* 56.23 272.39 1.001 3285.34 3390.05 1
Intercept[5] 4344.29* 1212.85 21046.56 1.001 3595.34 3536.92 1
Within-Person Effects
Daily persuasion experienced 0.84* 0.71 0.99 1.001 8449.12 8008.36 0.98
Daily persuasion utilized (partner’s view) 1.02 0.83 1.25 1.000 8067.47 7409.36 0.6
Daily pressure experienced 1.87* 1.17 2.77 1.001 5678.61 6449.73 0.994
Daily pressure utilized (partner’s view) 1.24 0.72 2.10 1.000 7142.04 6035.51 0.812
Daily pushing experienced 1.17 0.96 1.45 1.000 6912.78 7372.54 0.944
Daily pushing utilized (partner’s view) 0.91 0.70 1.18 1.000 8461.27 8004.37 0.76
Day 1.48 0.73 2.98 1.000 12423.41 8896.43 0.864
Daily weartime NA NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 1.14 0.40 3.33 1.000 4035.93 5849.76 0.594
Mean persuasion utilized (partner’s view) 1.43 0.46 4.69 1.000 4047.45 6116.21 0.728
Mean pressure experienced 3.56* 1.16 11.26 1.000 4687.35 6068.94 0.986
Mean pressure utilized (partner’s view) 1.16 0.35 3.63 1.000 4389.44 6592.04 0.6
Mean pushing experienced 1.21 0.27 5.75 1.001 4930.89 7590.39 0.593
Mean pushing utilized (partner’s view) 0.10* 0.01 0.64 1.001 6509.30 7214.09 0.995
Mean weartime NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.85 0.49 1.32 1.00 3877.68 5300.12 NA
sd(Daily persuasion experienced) 0.18 0.01 0.43 1.00 2038.71 4044.00 NA
sd(Daily persuasion utilized (partner’s view)) 0.22 0.01 0.52 1.00 3139.14 4498.36 NA
sd(Daily pressure experienced) 0.59 0.09 1.20 1.00 2757.85 2612.71 NA
sd(Daily pressure utilized (partner’s view)) 0.52 0.02 1.59 1.00 2807.86 4173.72 NA
sd(Daily pushing experienced) 0.22 0.01 0.53 1.00 3088.60 3994.39 NA
sd(Daily pushing utilized (partner’s view)) 0.19 0.01 0.62 1.00 4619.62 5640.25 NA
Additional Parameters
ar[1] 0.00 -0.92 0.94 1.00 3254.41 5617.26 NA
ma[1] NA NA NA NA NA NA NA
cosy NA NA NA NA NA NA NA
nu NA NA NA NA NA NA NA
shape NA NA NA NA NA NA NA
sderr 0.34 0.01 1.00 1.00 887.95 1136.18 NA
sigma NA NA NA NA NA NA NA
disc 1.00 1.00 1.00 NA NA NA NA 
hypothesis(reactance_ordinal, "pressure_self_cw > pushing_self_cw")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (pressure_self_cw... > 0     0.47      0.25     0.03     0.85      22.81
##   Post.Prob Star
## 1      0.96    *
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Binary

introduce_binary_reactance <- function(data) {
  data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
  return(data)
}



df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
  for (i in seq_along(implist)) {
    implist[[i]] <- introduce_binary_reactance(implist[[i]])
  }
}


formula <- bf(
  is_reactance ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw +
    
    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    pushing_self_cb + pushing_partner_cb +
    
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +
    pushing_self_cw + pushing_partner_cw | coupleID)
  , autocor = autocor_str
  )



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5) # range of the outcome scale
  
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  
  #, brms::set_prior("cauchy(0, 2)", class = "sigma", lb = 0)
  #, brms::set_prior("cauchy(0, 2)", class = "sderr", lb = 0)
  
  , autocor_prior
)


brms::validate_prior(
  prior1, 
  formula = formula, 
  data = df_double, 
  family = bernoulli()
  )
## Warning: Rows containing NAs were excluded from the model.
#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

is_reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::bernoulli(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("is_reactance_noar_", suffix))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
check_brms(is_reactance)
## 
## Divergences:
## 0 of 12000 iterations ended with a divergence.
## 
## Tree depth:
## 0 of 12000 iterations saturated the maximum tree depth of 10.
## 
## Energy:
## E-BFMI indicated no pathological behavior.

## Using 10 posterior draws for ppc type 'ecdf_overlay' by default.

## Using 10 posterior draws for ppc type 'dens_overlay' by default.

## Warning: Found 721 observations with a pareto_k > 0.7 in model 'model'. We
## recommend to set 'moment_match = TRUE' in order to perform moment matching for
## problematic observations.
## 
## Computed from 12000 by 756 log-likelihood matrix.
## 
##          Estimate   SE
## elpd_loo   -328.4 12.7
## p_loo       260.2 10.7
## looic       656.8 25.5
## ------
## MCSE of elpd_loo is NA.
## MCSE and ESS estimates assume MCMC draws (r_eff in [0.6, 1.2]).
## 
## Pareto k diagnostic values:
##                          Count Pct.    Min. ESS
## (-Inf, 0.7]   (good)      35    4.6%   396     
##    (0.7, 1]   (bad)      572   75.7%   <NA>    
##    (1, Inf)   (very bad) 149   19.7%   <NA>    
## See help('pareto-k-diagnostic') for details.
DHARMa.check_brms.all(is_reactance, integer = FALSE)

## Object of Class DHARMa with simulated residuals based on 1000 simulations with refit = FALSE . See ?DHARMa::simulateResiduals for help. 
##  
## Scaled residual values: 0.1160033 0.457393 0.4852882 0.4772632 0.05895521 0.9779663 0.0640845 0.7148414 0.9686449 0.5012651 0.6509825 0.2305564 0.9945411 0.01570549 0.687121 0.633468 0.6006051 0.4089159 0.4803867 0.1582965 ...
summarize_brms(
  is_reactance, 
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T,
  include_p_direction = T) %>%
  print_df(rows_to_pack = rows_to_pack)
OR l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS p_direction
Intercept 0.03* 0.00 0.18 1.001 2961.20 5776.43 1
Within-Person Effects
Daily persuasion experienced 0.55 0.26 1.01 1.000 4130.22 6033.87 0.972
Daily persuasion utilized (partner’s view) 1.86 0.76 6.03 1.001 3737.98 5220.40 0.908
Daily pressure experienced 8.70* 1.69 60.77 1.000 3710.48 6299.81 0.995
Daily pressure utilized (partner’s view) 2.19 0.33 18.42 1.000 6729.25 6651.90 0.789
Daily pushing experienced 2.34* 1.07 6.63 1.000 3633.87 4827.24 0.983
Daily pushing utilized (partner’s view) 0.72 0.21 2.27 1.001 6444.05 6815.31 0.723
Day 3.78 0.38 43.96 1.000 6637.82 7554.77 0.871
Daily weartime NA NA NA NA NA NA NA
Between-Person Effects
Mean persuasion experienced 13.42 0.90 291.21 1.000 3985.43 6196.30 0.97
Mean persuasion utilized (partner’s view) 5.05 0.31 91.34 1.000 7444.78 7562.53 0.872
Mean pressure experienced 132.06* 3.62 5363.94 1.000 9338.38 8352.40 0.997
Mean pressure utilized (partner’s view) 6.03 0.15 270.01 1.000 8879.00 9243.83 0.824
Mean pushing experienced 4.47 0.10 220.19 1.000 7718.10 8220.27 0.781
Mean pushing utilized (partner’s view) 0.06 0.00 3.21 1.000 8490.07 8668.77 0.92
Mean weartime NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 3.39 1.86 5.30 1.00 2132.54 3149.41 NA
sd(Daily persuasion experienced) 0.58 0.03 1.50 1.00 1978.49 4478.41 NA
sd(Daily persuasion utilized (partner’s view)) 1.50 0.42 2.95 1.00 2329.46 2696.99 NA
sd(Daily pressure experienced) 1.97 0.11 4.47 1.00 1972.68 3600.22 NA
sd(Daily pressure utilized (partner’s view)) 1.35 0.05 3.81 1.00 4817.26 5523.96 NA
sd(Daily pushing experienced) 0.84 0.05 2.02 1.00 2576.31 3993.09 NA
sd(Daily pushing utilized (partner’s view)) 0.75 0.03 2.23 1.00 4299.27 6263.29 NA
Additional Parameters
ar[1] 0.06 -0.16 0.27 1.00 1980.35 3211.48 NA
ma[1] NA NA NA NA NA NA NA
cosy NA NA NA NA NA NA NA
nu NA NA NA NA NA NA NA
shape NA NA NA NA NA NA NA
sderr 5.90 3.10 9.31 1.00 1442.76 1679.03 NA
sigma NA NA NA NA NA NA NA 
hypothesis(is_reactance, "pressure_self_cw > pushing_self_cw")
## Hypothesis Tests for class b:
##                 Hypothesis Estimate Est.Error CI.Lower CI.Upper Evid.Ratio
## 1 (pressure_self_cw... > 0     1.31      0.99    -0.23     3.01      11.11
##   Post.Prob Star
## 1      0.92     
## ---
## 'CI': 90%-CI for one-sided and 95%-CI for two-sided hypotheses.
## '*': For one-sided hypotheses, the posterior probability exceeds 95%;
## for two-sided hypotheses, the value tested against lies outside the 95%-CI.
## Posterior probabilities of point hypotheses assume equal prior probabilities.

Report All Models

if (report_ordinal) {
  model_rows_random_final <- model_rows_random_ordinal
  model_rows_fixed_final <- model_rows_fixed_ordinal
  model_rownames_fixed_final <- model_rownames_fixed_ordinal
  model_rownames_random_final <- model_rownames_random_ordinal
  rows_to_pack_final <- rows_to_pack_ordinal
} else {
  model_rows_random_final <- model_rows_random_hu
  model_rows_fixed_final <- model_rows_fixed_hu
  model_rownames_fixed_final <- model_rownames_fixed_hu
  model_rownames_random_final <- model_rownames_random_hu
  rows_to_pack_final <- rows_to_pack_hu
}



all_models <- report_side_by_side(
  pa_sub,
  pa_obj_log,
  mood,
  reactance,
  reactance_ordinal,
  is_reactance,
  
  model_rows_random = model_rows_random_final,
  model_rows_fixed = model_rows_fixed_final,
  model_rownames_random = model_rownames_random_final,
  model_rownames_fixed = model_rownames_fixed_final
) 
## [1] "pa_sub"
## [1] "pa_obj_log"
## Warning in summarize_brms(model, short_version = TRUE, exponentiate =
## exponentiate, : Coefficients were exponentiated. Double check if this was
## intended.
## [1] "mood"
## [1] "reactance"
## Warning: There were 330 divergent transitions after warmup. Increasing
## adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup
## Warning: Coefficients were exponentiated. Double check if this was intended.
## [1] "reactance_ordinal"
## [1] "is_reactance"
# pretty printing
summary_all_models <- all_models %>%
  print_df(rows_to_pack = rows_to_pack_final) %>%
  add_header_above(
    c(" ", "Subjective MVPA Hurdle Lognormal" = 2, 
      "Device-Based MVPA Lognormal" = 2, 
      "Mood Skewnormal" = 2,
      "Reactance Lognormal" = 2, 
      "Reactance Ordinal" = 2,
      "Reactance Dichotome" = 2)
  )

export_xlsx(
  summary_all_models, 
  rows_to_pack = rows_to_pack_final,
  file.path("Output", "AllModels_experimental.xlsx"), 
  merge_option = 'both', 
  simplify_2nd_row = TRUE,
  colwidths = c(38, 7.2, 13.3, 7.2, 13.3, 7.2, 13.3,7.2, 13.3,7.2, 13.3,7.2, 13.3),
  line_above_rows = c(1,2),
  line_below_rows = c(-1)
)
## 
## Attaching package: 'rvest'
## The following object is masked from 'package:readr':
## 
##     guess_encoding
summary_all_models
Subjective MVPA Hurdle Lognormal
Device-Based MVPA Lognormal
Mood Skewnormal
Reactance Lognormal
Reactance Ordinal
Reactance Dichotome
b pa_sub 95% CI pa_sub exp(Est.) pa_obj_log 95% CI pa_obj_log b mood 95% CI mood exp(Est.) reactance 95% CI reactance OR reactance_ordinal 95% CI reactance_ordinal OR is_reactance 95% CI is_reactance
Intercept 3.89* [ 3.76, 4.01] 117.92* [105.70, 131.58] 3.73* [ 3.58, 3.88] 1.80* [1.49, 2.15] NA NA 0.03* [0.00, 0.18]
Hurdle Intercept 0.18 [-0.15, 0.50] NA NA NA NA 4.59* [2.34, 9.06] NA NA NA NA
Conditional Within-Person Effects
Daily persuasion experienced 0.03 [-0.02, 0.08] 1.03 [ 1.00, 1.06] -0.01 [-0.03, 0.02] 0.99 [0.92, 1.07] 0.84* [0.71, 0.99] 0.55 [0.26, 1.01]
Daily persuasion utilized (partner’s view) 0.03 [-0.01, 0.08] 1.02 [ 0.99, 1.05] 0.01 [-0.01, 0.04] 1.00 [0.93, 1.06] 1.02 [0.83, 1.25] 1.86 [0.76, 6.03]
Daily pressure experienced -0.10* [-0.20, -0.01] 0.95 [ 0.89, 1.01] -0.02 [-0.09, 0.05] 1.06 [0.93, 1.18] 1.87* [1.17, 2.77] 8.70* [1.69, 60.77]
Daily pressure utilized (partner’s view) -0.07 [-0.16, 0.02] 0.98 [ 0.92, 1.04] -0.02 [-0.10, 0.05] 1.05 [0.87, 1.29] 1.24 [0.72, 2.10] 2.19 [0.33, 18.42]
Daily pushing experienced 0.03 [-0.04, 0.10] 1.03 [ 0.98, 1.07] 0.03 [-0.01, 0.07] 1.01 [0.95, 1.08] 1.17 [0.96, 1.45] 2.34* [1.07, 6.63]
Daily pushing utilized (partner’s view) -0.01 [-0.07, 0.06] 1.03 [ 0.99, 1.07] 0.04 [-0.01, 0.08] 1.01 [0.90, 1.15] 0.91 [0.70, 1.18] 0.72 [0.21, 2.27]
Day -0.02 [-0.15, 0.11] 0.96 [ 0.88, 1.05] 0.08 [ 0.00, 0.16] 1.00 [0.78, 1.28] 1.48 [0.73, 2.98] 3.78 [0.38, 43.96]
Daily weartime NA NA 1.00* [ 1.00, 1.00] NA NA NA NA NA NA NA NA
Conditional Between-Person Effects
Mean persuasion experienced 0.03 [-0.29, 0.35] 1.09 [ 0.82, 1.45] 0.23 [-0.17, 0.63] 0.94 [0.65, 1.36] 1.14 [0.40, 3.33] 13.42 [0.90, 291.21]
Mean persuasion utilized (partner’s view) 0.02 [-0.30, 0.33] 0.96 [ 0.71, 1.29] 0.06 [-0.34, 0.46] 1.22 [0.82, 1.85] 1.43 [0.46, 4.69] 5.05 [0.31, 91.34]
Mean pressure experienced 0.09 [-0.29, 0.45] 0.97 [ 0.71, 1.34] -0.13 [-0.52, 0.26] 1.21 [0.87, 1.72] 3.56* [1.16, 11.26] 132.06* [3.62, 5363.94]
Mean pressure utilized (partner’s view) -0.12 [-0.50, 0.27] 0.98 [ 0.72, 1.33] -0.20 [-0.58, 0.20] 0.86 [0.60, 1.23] 1.16 [0.35, 3.63] 6.03 [0.15, 270.01]
Mean pushing experienced 0.27 [-0.20, 0.75] 1.01 [ 0.66, 1.54] 0.02 [-0.54, 0.59] 0.85 [0.52, 1.36] 1.21 [0.27, 5.75] 4.47 [0.10, 220.19]
Mean pushing utilized (partner’s view) 0.27 [-0.21, 0.76] 1.29 [ 0.85, 1.96] 0.37 [-0.19, 0.94] 1.00 [0.47, 2.04] 0.10* [0.01, 0.64] 0.06 [0.00, 3.21]
Mean weartime NA NA 1.00 [ 1.00, 1.00] NA NA NA NA NA NA NA NA
Hurdle Within-Person Effects
Hu Daily persuasion experienced -0.43* [-0.56, -0.31] NA NA NA NA 1.19 [0.98, 1.47] NA NA NA NA
Hu Daily persuasion utilized (partner’s view) -0.28* [-0.40, -0.17] NA NA NA NA 0.96 [0.71, 1.31] NA NA NA NA
Hu Daily pressure experienced 0.21 [-0.12, 0.54] NA NA NA NA 0.53 [0.26, 1.06] NA NA NA NA
Hu Daily pressure utilized (partner’s view) -0.41* [-0.86, -0.05] NA NA NA NA 0.75 [0.25, 2.25] NA NA NA NA
Hu Daily pushing experienced -0.55* [-0.92, -0.24] NA NA NA NA 0.80 [0.63, 1.01] NA NA NA NA
Hu Daily pushing utilized (partner’s view) -0.62* [-0.89, -0.39] NA NA NA NA 1.14 [0.78, 1.75] NA NA NA NA
Hu Day 0.09 [-0.17, 0.35] NA NA NA NA 0.61 [0.28, 1.30] NA NA NA NA
Hu Daily weartime NA NA NA NA NA NA NA NA NA NA NA NA
Hurdle Between-Person Effects
Hu Mean persuasion experienced -0.19 [-1.05, 0.67] NA NA NA NA 0.63 [0.15, 2.32] NA NA NA NA
Hu Mean persuasion utilized (partner’s view) -0.18 [-1.02, 0.67] NA NA NA NA 0.59 [0.13, 2.53] NA NA NA NA
Hu Mean pressure experienced 1.25* [ 0.33, 2.16] NA NA NA NA 0.05* [0.00, 0.50] NA NA NA NA
Hu Mean pressure utilized (partner’s view) 0.61 [-0.30, 1.51] NA NA NA NA 0.84 [0.08, 10.66] NA NA NA NA
Hu Mean pushing experienced -1.14 [-2.36, 0.09] NA NA NA NA 1.75 [0.21, 17.49] NA NA NA NA
Hu Mean pushing utilized (partner’s view) -1.16 [-2.39, 0.08] NA NA NA NA 25.69* [2.54, 336.83] NA NA NA NA
Hu Mean weartime NA NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.31 [0.24, 0.41] 0.30 [0.23, 0.39] 0.42 [0.33, 0.55] 0.17 [ 0.01, 0.36] 0.85 [ 0.49, 1.32] 3.39 [ 1.86, 5.30]
sd(Hurdle Intercept) 0.90 [0.70, 1.17] NA NA NA NA 1.12 [ 0.67, 1.73] NA NA NA NA
sd(Daily persuasion experienced) 0.11 [0.07, 0.16] 0.05 [0.03, 0.08] 0.02 [0.00, 0.05] 0.11 [ 0.04, 0.20] 0.18 [ 0.01, 0.43] 0.58 [ 0.03, 1.50]
sd(Daily persuasion utilized (partner’s view)) 0.08 [0.04, 0.13] 0.05 [0.02, 0.08] 0.02 [0.00, 0.05] 0.06 [ 0.00, 0.16] 0.22 [ 0.01, 0.52] 1.50 [ 0.42, 2.95]
sd(Daily pressure experienced) 0.06 [0.00, 0.20] 0.04 [0.00, 0.13] 0.04 [0.00, 0.14] 0.12 [ 0.01, 0.28] 0.59 [ 0.09, 1.20] 1.97 [ 0.11, 4.47]
sd(Daily pressure utilized (partner’s view)) 0.06 [0.00, 0.18] 0.04 [0.00, 0.11] 0.06 [0.00, 0.17] 0.15 [ 0.01, 0.47] 0.52 [ 0.02, 1.59] 1.35 [ 0.05, 3.81]
sd(Daily pushing experienced) 0.12 [0.06, 0.20] 0.06 [0.00, 0.14] 0.03 [0.00, 0.07] 0.07 [ 0.00, 0.19] 0.22 [ 0.01, 0.53] 0.84 [ 0.05, 2.02]
sd(Daily pushing utilized (partner’s view)) 0.11 [0.04, 0.18] 0.03 [0.00, 0.08] 0.04 [0.00, 0.10] 0.12 [ 0.01, 0.29] 0.19 [ 0.01, 0.62] 0.75 [ 0.03, 2.23]
sd(Hu Daily persuasion experienced) 0.18 [0.02, 0.34] NA NA NA NA 0.22 [ 0.01, 0.51] NA NA NA NA
sd(Hu Daily persuasion utilized (partner’s view)) 0.17 [0.02, 0.33] NA NA NA NA 0.49 [ 0.08, 0.99] NA NA NA NA
sd(Hu Daily pressure experienced) 0.30 [0.01, 0.85] NA NA NA NA 1.11 [ 0.13, 2.45] NA NA NA NA
sd(Hu Daily pressure utilized (partner’s view)) 0.35 [0.01, 1.02] NA NA NA NA 1.10 [ 0.05, 3.22] NA NA NA NA
sd(Hu Daily pushing experienced) 0.65 [0.33, 1.08] NA NA NA NA 0.24 [ 0.01, 0.59] NA NA NA NA
sd(Hu Daily pushing utilized (partner’s view)) 0.32 [0.05, 0.65] NA NA NA NA 0.32 [ 0.01, 0.99] NA NA NA NA
Additional Parameters
ar[1] 0.35 [0.27, 0.42] 0.29 [0.26, 0.33] 0.90 [0.74, 1.00] 0.19 [-0.85, 0.94] 0.00 [-0.92, 0.94] 0.06 [-0.16, 0.27]
ma[1] NA NA NA NA NA NA NA NA NA NA NA NA
cosy NA NA NA NA NA NA NA NA NA NA NA NA
nu NA NA NA NA NA NA NA NA NA NA NA NA
shape NA NA NA NA NA NA NA NA NA NA NA NA
sderr 0.62 [0.60, 0.65] NA NA 0.20 [0.15, 0.25] 0.19 [ 0.01, 0.38] 0.34 [ 0.01, 1.00] 5.90 [ 3.10, 9.31]
sigma NA NA 0.55 [0.54, 0.57] 0.86 [0.82, 0.90] 0.33 [ 0.18, 0.44] NA NA NA NA 
report::report_system()

Analyses were conducted using the R Statistical language (version 4.4.1; R Core Team, 2024) on Windows 11 x64 (build 22635)

report::report_packages()
  • beepr (version 2.0; Bååth R, 2024)
  • R.methodsS3 (version 1.8.2; Bengtsson H, 2003)
  • R.oo (version 1.26.0; Bengtsson H, 2003)
  • R.utils (version 2.12.3; Bengtsson H, 2023)
  • brms (version 2.21.0; Bürkner P, 2017)
  • Rcpp (version 1.0.13; Eddelbuettel D et al., 2024)
  • bayesplot (version 1.11.1; Gabry J, Mahr T, 2024)
  • lubridate (version 1.9.3; Grolemund G, Wickham H, 2011)
  • DHARMa (version 0.4.6; Hartig F, 2022)
  • wbCorr (version 0.1.22; Küng P, 2023)
  • tibble (version 3.2.1; Müller K, Wickham H, 2023)
  • R (version 4.4.1; R Core Team, 2024)
  • openxlsx (version 4.2.7.1; Schauberger P, Walker A, 2024)
  • ggplot2 (version 3.5.1; Wickham H, 2016)
  • forcats (version 1.0.0; Wickham H, 2023)
  • stringr (version 1.5.1; Wickham H, 2023)
  • rvest (version 1.0.4; Wickham H, 2024)
  • tidyverse (version 2.0.0; Wickham H et al., 2019)
  • readxl (version 1.4.3; Wickham H, Bryan J, 2023)
  • dplyr (version 1.1.4; Wickham H et al., 2023)
  • purrr (version 1.0.2; Wickham H, Henry L, 2023)
  • readr (version 2.1.5; Wickham H et al., 2024)
  • xml2 (version 1.3.6; Wickham H et al., 2023)
  • tidyr (version 1.3.1; Wickham H et al., 2024)
  • knitr (version 1.48; Xie Y, 2024)
  • kableExtra (version 1.4.0; Zhu H, 2024)
report::cite_packages()